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Experiments on isotropic compression of a granular assembly of spheres show that the shear and 
bulk moduli vary with the confining pressure faster than the 1/3 power law predicted by Hertz- 
Mindlin effective medium theories (EMT) of contact elasticity. Moreover, the ratio between the 
moduli is found to be larger than the prediction of the elastic theory by a constant value. The 
understanding of these discrepancies has been a longstanding question in the field of granular matter. 
Here we perform a test of the applicability of elasticity theory to granular materials. We perform 
sound propagation experiments, numerical simulations and theoretical studies to understand the 
elastic response of a deforming granular assembly of soft spheres under isotropic loading. Our 
results for the behavior of the elastic moduli of the system agree very well with experiments. We 
show that the elasticity partially describes the experimental and numerical results for a system under 
compressional loads. However, it drastically fails for systems under shear perturbations, particularly 
for packings without tangential forces and friction. Our work indicates that a correct treatment 
should include not only the purely elastic response but also collective relaxation mechanisms related 
to structural disorder and nonaffine motion of grains. 



I. INTRODUCTION AND OBJECTIVES 

The acoustic properties of granular porous materials 
confined by an external stress can be extremely nonlin- 
ear as compared to continuum elastic solids 0, 0, 0, • 
Many industrial applications, such as the optimization 
of well location in an oil reservoir, depend crucially on 
the correct interpretation of nonlinear acoustic effects in 
granular materials, as exemplified by the large variation 
of the sound speeds or the elastic constants of the gran- 
ular formation as a function of the external stress |jj E3 ■ 

Important insight into this problem comes first from 
the Hertz-Mindlin contact theory to model the intergrain 
forces 0, 13- In this case, nonlinearity arises from the 
increase with the external stress of the contact area be- 
tween two spherical grains. Conventional theories de- 
scribing this problem in the framework of elasticity of 
continuum media |8( consider a uniform strain at all 
scales, and the displacement field of the grains is affinc 
with the macroscopic deformation (the affine approxima- 
tion). Here, one computes the stresses in terms of the 
strains by considering the disordered medium as an ef- 
fective medium that exerts a mean-field force (as given 
by contact Hertzian theory) on a single representative 
grain. This approximation is usuall y re ferred to as the 
Effective Medium Theory (EMT) [IM El 113 

As shown in a short letter E3 and the studies of other 
groups , the EMT does not successfully explain the me- 
chanical properties of cohesionless granular assemblies. 
The main prediction of the theory is the scaling of the 
bulk modulus, K, and shear modulus, /i, with the pres- 
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However there is a large 
volume of experiments for irregular sand grains as well 
as spherical glass beads which show anomalous scaling 
characterized by exponents varying between 1 /3 and 1 /2 



(for a comprehensive review see Goddard [3 and for a re- 
view in the geotechnical literature see Some studies 
have suggested that a ~ p 1 ^ 2 scaling is more appropriate 
for describing the nonlinear variation of the moduli . 

Here we extend the results of E3 an d investigate the 
applicability of elasticity theory to granular matter by 
means of experiments, computer simulations, and ana- 
lytical calculations. We first develop a series of acous- 
tic experiments to characterize the nonlinear elastic be- 
havior of non-cohesive dry granular materials under a 
wide range of external pressures. From this experimen- 
tal study we conclude that a microscopic study is needed 
in order to elucidate the deficiencies of existing granu- 
lar theories. Then we perform a Molecular Dynamics 
(MD) simulation to give microscopic insight into the re- 
laxation mechanism of granular materials. Finally, we 
offer a simple theory of relaxation going beyond the ef- 
fective medium approximation of elasticity. 

We calculate the elastic moduli, K and /i, of a dis- 
ordered array of elasto-frictional Hertz-Mindlin spherical 
grains. Numerical simulations resolve the question as to 
whether the problem lies in the treatment of intergrain 
contact or with the EMT. We find agreement between our 
simulations and the experiments, thus confirming the va- 
lidity of the Hertz-Mindlin contact theory to glass bead 
aggregates composed of frictional particles. 

Regarding the anomalous pressure dependence of the 
moduli, we find that there are several nonlinearities 
which preclude the proper definition of a scaling behav- 
ior as a function of pressure. We find a regime at low 
pressure where the coordination number and the volume 
fraction do not change much from their minimal values. 
In this regime the p 1 ^ 3 scaling is approximately valid. 
However, for pressures larger than 10 MPa the increase 
of the coordination number and volume fraction induces 
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other nonlinearities and therefore no simple scaling be- 
havior can be defined. 

In reality around 10 MPa there is crossover from a p 1 / 3 
to a p 5 / 9 scaling at larger pressures. Thus, we conclude 
that the scattering in the experimental values of the pres- 
sure exponent might be explained by the fact that the 
exponent is actually changing continuously from 1/3 to 
5/9. This is especially true in the regime where the ex- 
periments are usually done, near the crossover region at 
10 MPa. Similar conclusions have been reached by Lud- 
ing who also found that EMT needs to take into 
account the coordination number and its dependence on 
pressure. 

Our most important results relate to the effect of fric- 
tion and stress relaxation on the behavior of the elas- 
tic constants. Firstly, we find that the elastic formula- 
tion gives a reasonable description for the response of the 
system to compressional loads, i.e. the bulk modulus is 
reasonably well defined with the simple EMT. However, 
our simulations establish that the EMT is inadequate in 
describing the response of the material under shear per- 
turbations. 

The numerical simulations indicate that EMT fails be- 
cause it does not properly allow the grains to relax from 
the affine deformation imposed by the external bound- 
ary. The affine assumption is that, under an infinitesimal 
symmetric macroscopic deformation, each grain trans- 
lates according to the direction of the macroscopic strain, 
and it does not rotate. Moreover, no further relaxation 
mechanism is allowed. Such an homogeneous strain field 
is consistent with the local force balance of grains only 
in an ordered system. For disordered packings an inho- 
mogeneous strain develops at the local level. After the 
application of an affine strain the particles experience 
an unbalanced force since they are not, in general, in a 
symmetric environment. Consequently, the particles will 
move to a position different to that predicted by the affine 
approximation, so that the net force on each particle is 
zero. Similarly for torques and rotations. 

Here we show that the assumption of affinity is approx- 
imately valid for the bulk modulus and seriously flawed 
for the shear modulus. For this reason, the EMT pre- 
diction differs significantly from the experimental value. 
Thus, the principal source of deviation from EMT is the 
breakdown of the uniform strain assumption. 

To quantify the breakdown of EMT for the shear mod- 
ulus we focus our studies to two cases: frictionless grains 
interacting via only normal forces (this system is said to 
be path-independent and it is thought to describe com- 
pressible foams and emulsions) and systems with elas- 
tic tangential forces and Coulomb frictional forces (these 
systems are path-dependent and describe dry granular 
materials). 

The largest disagreement between theory and simula- 
tions is found for frictionless systems; the difference is 
more pronounced at low confining stress where we show 
that the system is in a state of marginal rigidity at a 
minimal mean coordination number equal to 6 in 3D (or 



4 in 2D). We find that after the application of an ex- 
ternal shear strain there is a nearly complete relaxation 
of the system to the applied shear; a result that cannot 
be captured under the framework of elasticity. We show 
that a new scaling behavior with pressure might describe 
the data for frictionless particles better: p(p) ~ p 2 / 3 as 
p -> 0. 

We interpret this result in the framework of critical 
phenomena: as p — > the system approaches a critical 
point at a mean coordination number Z c = 2D in D di- 
mensions, and a volume fraction of random close packing 
4> c w 0.64. This point describes a rigidity threshold state 
or a critical state of the packing as defined by Alexan- 
der ml and it is where the system becomes "isostatic" 
[Til Il8l Il9| . The elastic moduli vanish as a power-law 
of the pressure or volume fraction. For any finite pres- 
sure rigidity is achieved, since Z > Z c . Near the rigidity 
threshold the reference packing structure has a power law 
dependence on the pressure, modifying the scaling of fi(p) 
predicted by the Hertz theory. 

When friction and tangential elasticity is restored at 
the intergrain contacts, the agreement between theory 
and simulations (and in this case experiments) improves 
with respect to the frictionless case. This is because the 
existence of tangential restoring forces reduces the ex- 
tent to which the grains relax from the assumed affine 
configuration. Thus, the EMT provides a better agree- 
ment with simulations and with experiments for frictional 
grains than for frictionless grains, but serious disagree- 
ments still persist as we shall demonstrate. 

We conclude that in order to develop a better under- 
standing of the problem, one must abandon the purely 
elastic framework and consider granular matter as a full 
viscoelastic body. Collective relaxation effects can ac- 
count for the discrepancy in the shear modulus in com- 
parison with the elastic prediction: the corrections in- 
crease dramatically in the case of loose materials and for 
frictionless packings. A theory of single-particle relax- 
ation is offered as a first step in this direction. We also 
discuss our results in the framework of recent theories of 
marginal rigidity, jamming, melting and fragile matter. 

Applications. - Part of the motivation for this research 
derives from the fact that acoustics and nonlinear elas- 
tic logging methods are at the forefront of the evolving 
technology to help plan and optimize well location in oil 
exploration |20j. In order to position a well correctly, 
the knowledge of the stress distribution around the bore- 
hole is essential. Mechanical properties of the granular 
formation obtained from sonic logging can help predict 
formation strength, while stress magnitude derived from 
sonic measurements helps in predicting sanding problems 
in unconsolidated formations. Acoustic measurements in 
granular materials provide the natural way to understand 
the distribution of stress around the borehole. 

Quite apart from the relevance to borehole logging, 
within the field of seismic tomography there is a grow- 
ing interest in developing techniques to generate images 
of dissipation, along with the more traditional images of 
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impedance contrast. This extended seismic tomography 
would have impact, not only on the understanding of hy- 
drocarbon reservoirs but also e.g. on techniques to mon- 
itor the migration of ground water pollutants. This work 
represents an attempt to understand some of the mecha- 
nisms of attenuation as well as the stress dependence of 
sound velocities in granular materials. 

The outline of the paper is as follows. Section [H] re- 
views the background of the problem of Effective Medium 
theories and numerical approaches: MD simulations and 
intergrain forces for granular materials, compressed emul- 
sions and foams. Section IlIII describes the experiments 
on sound propagation. Section llVI describes the numer- 
ical results and Section [V] t ne theoretical results. We 
conclude in Section IvU with a final outlook. 



= when £ < 0. The effective 
Vg) is defined in terms of the 



II. BACKGROUND 

The problem of elastic properties of granular materials 
has been treated by many researchers since the pio nee ring 
work of Mindlin in the 50 s 0, 0, H3, [53, [53, HE Ell 
[2^, [2j|[3(j. However, a general solution to this problem 
is still lacking. 

In a typical experiment, a set of cohesionless glass 
beads is confined at a hydrostatic stress, p, and the com- 
pressional sound speed, v p , and the shear sound speed, 
v s , are measured as functions of stress (see for instance 
Domenico H3, Yin [H, and 0, [H H|). The P-wave 
and S-wave speeds are related to the elastic constants of 
the material in the long- wavelength limit: 



'K + 4/3/J, 



(1) 



(2) 



where p is the mass density of the system. 



A. Contact mechanics. 

In his seminal paper "On the Contact of Elastic Solids" 
H. Hertz 0, @ used linear elasticity of continuum me- 
dia to calculate the normal force of two perfectly elas- 
tic spheres pressed into contact considering no attraction 
or stickiness. Hertz showed that two spherical grains in 
contact with radii R\ and R2 interact with a normal re- 
pulsive force 



only in compression, F n 
stiffness k n — 4fi g /(l - 
shear modulus of the grains p g and the Poisson ratio v g 
of the material from which the grains are made (typically 
H g = 29 GPa and v g = 0.2, for spherical glass beads). 

The situation in the presence of a tangential force, F t , 
is more complicated. In the case of spheres under oblique 
loading, the tangential contact force was first calculated 
by Mindlin [3lJ. A general loading history can be de- 
scribed by the incremental change in the tangential force 
AF t and in the normal force AF n . For the special case 
where the partial increments do not involve microslip at 
the contact surface (i.e., |AF t | < fifAF n , where pj is the 
kinematic friction coefficient between the spheres, typi- 
cally /if = 0.3) Mindlin [3lJ showed that the tangential 
force is 



AF t = k t (R0 1/2 As, 



(4) 



where k t = 8/j, s /(2 — v g ), and the variable s is defined 
such that the relative shear displacement between the two 
grain centers is 2s. This is called the Mindlin "no-slip" 
solution. [See Appendix A for a more general solution] . 

The incremental form Eq. (0} is needed since the nu- 
merical value of the tangential force depends upon the 
trajectory taken in the space (£, s), see for details. 
The tangential force is obtained by integrating over the 
path taken by the spheres in contact subject to the initial 
conditions: F n = 0, i*t=0at£ = 0,s = 0, yielding: 



Ft = 



k t (R£) 1/2 ds. 



(5) 



path 



Thus, a granular system with tangential elastic forces 
is said to be path-dependent. By path dependency we 
mean that the work done in deforming the system de- 
pends upon whether one first compresses the system, 
then shears it, or first shears it then compresses. The 
results depend upon the path taken and not just the in- 
stantaneous final state. On the other hand, a system 
of spheres interacting only via normal forces, Eq. |J3J|, 
is said to be path-independent, and the work does not 
depend on the way the strain is applied. 

As the shear displacement increases, the elastic tan- 
gential force F t reaches its limiting value given by Amon- 
tons' law for no adhesion, F t < /J,fF n . Amontons' law (a 
special case of Coulomb's law) adds a second source of 
path-dependency as well as hysteresis to the problem. 



Effective medium theories (EMT) of granular 
elasticity 



F n = ~ k n R^e /2 , (3) 

where R = 2R1R2KR1 + R2), the normal overlap is £ = 
(1/2) + -R2) — l^i — a?2 1] > 0, and X\, x 2 are the 
positions of the grain centers. The normal force acts 



The basic idea of elastic theories relevant to our study 
is that the macroscopic work done in deforming the sys- 
tem is set equal to the sum of the work done on each 
grain-grain contact and that the latter is replaced by a 
suitable average 0, E3- These theories are usu- 
ally referred to as the effective medium theory (EMT) 
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and are based on Hertz-Mindlin contact mechanics. In 
the case of an isotropic deformable solid (for simplicity 
we describe the isotropic case) , the strain energy density 
per unit volume as a function of the strain ey is: 



(6) 

This equation is equivalent to the expression 

dij = K e euSij + 2fi e (eij - -<%£«), (7) 

which determines the stress tensor cry in terms of the 
strain tensor for an isotropic body Here 

= l/2(dui/dxj + duj/dxi), (8) 

and the deviations Ui = X{ — Ri of the positions of the 
N particles in the system, {xi, . . . , xjv} , are measured 
from a suitable rigid reference state 



{R} — {Ri,. . . , Rat} , 



(9) 



around which one can expand consistently. This refer- 
ence state is straightforward for simple periodic systems. 
However, it is assumed that this expansion is also possi- 
ble for amorphous solids with reference states which are 
random, like in a packing of grains (see Alexander [l6| 
for extensive discussions). The subscript in K e and /i e 
denotes that the values of the moduli are calculated con- 
sidering granular materials as purely elastic solids. There 
are two assumptions inherent to the elastic EMT: 

• All the spheres are statistically the same, and it is 
assumed that there is an isotropic distribution of 
contacts around a given sphere. 

• An affine approximation is used, i.e., the spheres 
at position Xj are moved a distance 5ui in a time 
interval St according to the macroscopic strain rate 



by 



(10) 



The grains are always at equilibrium due to the 
assumption of an isotropic distribution of contacts, 
and further relaxation is not required. This sort of 
mean field theory is analogous to a simple average 
of non-linear spring constants. 

It is important to notice that a definition of uniform 
strain field is possible only under the mean field approx- 
imation. This assumption is also trivially correct for or- 
dered packings. However, for disordered systems, the 
affine approximation is inconsistent with the local equi- 
librium of grains . We will come back to this crucial 
point later on. 

With the above assumptions, the elastic energy Eq. 
© is set equal to a suitable average over the contacts, 
viz: 



U 



- E 

V ^ 

contact 



jF-du^^(jF-du), (11) 



with F du = F n d£,+F t -ds, Z is the average coordination 
number defined as the average number of contacts per 
particle, <fi is the volume fraction of the sample, and Vb 
is the volume of a single grain. The EMT predictions 
for the bulk and shear modulus for an isotropic system 
compressed at pressure p are the following: 
We distinguish between two different models: 
(1) Path-independent models, k t = 0, frictionless 
grains: We consider that there is perfect slippage at the 
intergrain contact. This corresponds to F t = 0, only nor- 
mal forces between the particles. This case corresponds 
to path-independent forces, and allows the use of an en- 
ergy density function Eq. © which depends only on 
the instantaneous position of the particles. This case 
could be considered conservative since the total work on 
a closed path is zero. A system of frictionless spherical 
particles could be thought of as a model of compressed 
emulsions and foams which are usually modeled as vis- 
coelastic spheres without tangential forces [3^, |3 HH 
(see Appendix C). 
For the case of frictionless grains one finds: 
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(12) 



(13) 



(2) Path- dependent models, kt 7^ 0, frictional grains: 
In this case tangential elastic forces are taken into con- 
sideration. In principle the energy functional © now de- 
pends on the path. However, it has been shown that the 
second order elastic constants are still path-independent 
under the framework of EMT, while path-dependency 
appears only in the third order elastic constants [l2|. 

The bulk modulus is not affected by the introduction 
of tangential forces, and Eq. (|12fl is still valid in this 
case. However, the shear modulus is modified according 
to: 



Me(p) 



k n 



20tt 
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(14) 



The above results have been obtained by a number 
of authors using different methods and are valid in 3-D 
[ToL ITU Il2| . It should be noted that the above results 
are obtained for a system of infinitely rough spheres, i.e. 
when fJ>f —> 00. Thus, there is no sliding, and Coulomb 
friction is not considered, although there is a tangential 
elastic restoring force as given by Eq. {3J. See Appendix 
A for further discussion. 

The p 1 / 3 dependence in Eqs. I|12|) - (|14|l is a direct con- 
sequence of the scaling of the normal Hertz force on the 
deformation. Since 



P 



3/2 3/2 



(15) 
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Then we expect 

Ve-K^^-e^-p 1 / 3 . (16) 

We note that a system of linear springs, F n ~ £, would 
give rise to elastic constants which are independent of 
pressure, as in the linear elasticity theory. 



that EMT does not correctly describe the shear modulus 
but it describes the bulk modulus fairly well. Other ex- 
perimental work done by Liu and Nagel 25] and Jia et 
al. [2(j concentrated on the role played by force chains 
in sound propagation. Different approaches for one di- 
mensional elastic chains have also been applied to wave 
propagation in granular media p9l l30l | . 



C. Discrepancies between theory and experiments. 

It was found experimentally that the shear and bulk 
moduli of an assembly of spherical grains vary with the 
confining stress, p, faster than the p 1 / 3 power law pre- 
dicted by Eqs. JTSJ-JTSJ 2]. Another way of seeing the 
breakdown of the elastic theory is to focus on the ra- 
tio K/fi. According to Eqs. I|12fl and (|14|) for frictional 
grains 

= 5(2 - v g ) ^ 
He 3(5 - Av g ) ' 

independent of stress, a value which depends only on the 
Poisson ratio of the bead material. The experiments give 
K/n w 1.1 - 1.3 0i|. EMT predicts K e /(i e = 0.71, 
if we take v g = 0.2 for the Poisson ratio of glass. [The 
EMT prediction is rather insensitive to variations of v g ; 
K e /ne = 0.71 ± 0.04 for v g = 0.2 ±0.1.] Conversely, a 
value v g ~ 1.2 would be needed in order to fit the experi- 
mental K/fi, clearly violating the upper thermodynamic 
limit of v g < 1/2 ||. 

Another quantity of interest is the effective Poisson ra- 
tio of the pack v. According to EMT, v e is again indepen- 
dent of pressure and given only in terms of the Poisson's 
ratio of the grains: 

def K e - 2/3He = Vg , . 

2(K e -l/3ne) 2(5-3^)" { ' 

Thus, for typical glass beads (v g — 0.2) we find a pre- 
dicted value of v e = 0.02 which is one order of magnitude 
smaller than typical experimental values v w 0.28 |2lj : 
another serious disagreement. 

The origin of the above discrepancies has not been 
clear: it could be due to the breakdown of the Hertz- 
Mindlin force law at each grain contact, or it could be 
associated with the breakdown of the elasticity theory 
applied to granular systems. De Gennes [27J proposed 
that a thin shell of oxide layer would give a faster growth 
with stress of the elastic moduli of the system, which 
may explain the behavior for metallic beads. Goddard 
proposed that sharp angularities of the grains (for in- 
stance sand grains) may modify the contact force law be- 
tween grains, giving rise to a different stress dependence. 
Other authors 0, have suggested that the increas- 
ing number of contacts with stress may be the reason for 
the discrepancies in the stress dependence of the moduli. 
Jenkins et al. |28j measured the elastic moduli using nu- 
merical simulations for a single pressure and concluded 



D. Linear viscoelastic constitutive models 

In order to understand our results it is important to 
generalize the elastic concepts introduced above to the 
full viscoelastic response. In linear viscoelasticity |36j . 
the current state of stress specified by the stress tensor cr^- 
is determined by the past history via a linear constitutive 
equation 

*ij (t) = I G m {t-t 1 ) e kl (O < (19) 

where e^i = deu/dt is the strain rate, and Gijki(t) is 
called the relaxation modulus tensor. 

For an isotropic linear viscoelastic material, the relax- 
ation modulus tensor has only two independent compo- 
nents. These are the shear relaxation modulus G(t) and 
the bulk relaxation modulus K(t) characterizing the re- 
sponse to shear £12 and bulk deformation in. The re- 
laxation modulus G{t) and K(t) are conceptualized as 
the time-dependent analogues of the shear \x e and bulk 
modulus K e in elasticity theory. 

In this study we will concentrate on the stress relax- 
ation after a sudden strain imposed via a simple shear, 
a pure shear, or a uniaxial compression. For example, a 
shear strain is applied instantaneously, at t = 0, from its 
initial value of zero to a final, constant value, £12- For 
this situation we have £12 (i) = ei2<K*)- Equation 119|) 
reduces to 

a 12 (t) =ei 2 G(t) . (20) 

Therefore, this strain protocol immediately yields com- 
plete information on the response function, G(t), [short- 
hand for Gi2i2(£) here] simply by measuring 012 (i)- This 
is a strain protocol which is particularly simple to imple- 
ment in our MD simulations. 

For a perfectly elastic solid the relaxation modulus is 
independent of time, G(t) = G ^constant, and one can 
define the shear modulus of the solid as 

Me = er 12 /ei2 = G. (21) 

We will show that the instantaneous response of the vis- 
coelastic granular material, G(t = 0), represents the 
shear modulus, /i e , as calculated by the effective medium 
theories of continuum elasticity. 

For a Newtonian liquid G(t) — 7]S(t), where r\ is the 
viscosity. For a viscoelastic liquid, G(t) approaches zero 
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as t — > oo. For a viscoelastic solid, structural relaxation 
and elasticity lead to a finite modulus as t — > oo: 

f i = G(t~>oo). (22) 

A similar analysis can be performed for the bulk modulus, 
K(t) denned as 

a ii (t) = 3e ii K(t) . (23) 
to obtain K(Q) = K e and 

K = K(t -> oo). (24) 

Equations (|22|l and (|24|) will be used to calculate the 
moduli in the simulations. 



E. Molecular Dynamics simulations 

In MD simulations of granular matter the net force 
and moment on each grain depend on the choice of in- 
tergrain contact laws |37l l3S||. Here, we follow the Dis- 
crete Element Method (DEM) developed by Cundall and 
Strack [37j and solve Newton's equations for an assembly 
composed of soft elasto-frictional spheres interacting via 
Hertz-Mindlin contact forces and Coulomb friction as de- 
scribed in Section III Al ■ We employ a time-stepping, 
finite-difference approach to solve the Newtonian equa- 
tions of motion simultaneously for every grain in the sys- 
tem: 

F = to x, (25) 



M = I 0, (26) 

where F and M are the net force and moment acting on a 
given grain, to and I are the mass and moment of inertia, 
and x and 9 are the linear and angular accelerations of 
the grain, respectively. 

The numerical solution of Eqs. (|25|l and (|26fl are ob- 
tained by integration, assuming constant velocities and 
accelerations for a given time step: linear and angular ve- 
locities are determined from the knowledge of the force 
and torque, and grain displacements and rotations at the 
next time step are calculated from the average velocities. 
Grain motions can be initiated by gravitational forces, by 
external forces prescribed by stress or strain rate bound- 
ary conditions, and by forces resolved at intergrain con- 
tacts. Strain rates are assumed to be low, and small time 
steps At are chosen to ensure that the disturbance of a 
given grain only propagates to its immediate neighbors 
(see Appendix D). 

Viscous damping. — Damping of grain motions must 
be included in the calculations to prevent the continuous 
oscillation of an elastic system. Although damping is a 
physical reality, and physically meaningful mechanisms 



might well be incorporated, our concern here is to get 
the simulations to equilibrate to the final answer in a 
reasonable amount of computer time. 

Several damping methods are possible. Global damp- 
ing considers the particles immersed in a viscous fluid 
and is provided by introducing viscous force terms in Eqs. 
(|25|l and l|26|) . These drag forces are proportional to the 
absolute velocity and angular velocity of the particles: 
~ — 7 n x, and ~ —7*0, where the 7's are the global damp- 
ing coefficient related to the viscosity of the immersing 
fluid (which could be for instance air). 

Global damping is introduced to guarantee that the 
system can reach an equilibrium state with zero veloc- 
ity at a given pressure. Its physical significance is be- 
ing studies at the moment by experiments and computer 
simulations. Another source of damping implies a con- 
tact force term acting at every contact point, propor- 
tional to the relative velocities of the grains. Microscopic 
contact damping occurs due to the viscous dissipation 
of energy in the bulk of the particle material when they 
are deformed and it may also occur if liquid bridges are 
formed at the contact points between the particles. Here, 
a damping force is added to each contact force, Eqs. © 
and J5J , proportional to the relative normal and shear ve- 
locities, /3 ra £ and fas, respectively, with (3 n and (3 t the con- 
tact damping coefficients. Typical values of the damping 
constant are given in |38l ]. 

In this study we will use global damping for the prepa- 
ration of the sample and the calculation of the elastic 
constants. This procedure is necessary to achieve the 
final equilibrium states which we wish to explore (see 
Appendix B for a discussion). 

Computation of stress. — The macroscopic stress ten- 
sor for point contacts in a volume V is given by [Tol ITU 

m 

a v = W 2 (FiRrij + RniFj), (27) 

contacts 

where n is the unit vector joining the center of two 
spheres in contact. 

III. ACOUSTIC EXPERIMENTS 

In the simplest experiments, a packing of glass beads is 
confined under hydrostatic conditions and the compres- 
sional and shear sound speeds, v p and v s , are measured 
as functions of p H El |H H3 . 

In the long- wavelength limit, the sound speeds are re- 
lated to the elastic constants of the aggregate by Eqs. (JTJ 
and J2J) • Here we perform our own experiments according 
to standard sound propagation techniques [2ll [2j| . 

A. Experimental Configuration 

We used a set of high quality glass beads of a small 
enough diameter to measure an appreciable signal at low 
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FIG. 1: Container and Transducers-LVDT apparatus used in 
the sound propagation experiments. 



FIG. 2: (a) Wave velocities versus pressure obtained in our 
experiments. Also shown are the results of Domenico |2l| 
for comparison. We cycle up and down in pressure to avoid 
hysteresis. 



pressure. From the experimental data of Domenico [2l| , 
we expect compressional velocities v p ~ 1000 m/s and 
shear velocities v s ~ 500 m/s at low pressures. We per- 
form ultrasonic measurements with pulses of frequency 
/ = 500 kHz, and we find that the maximum size of the 
beads should be R <C v s /(2f). Then, we choose a set 
of glass beads of diameter 45 /im in order to reach the 
desired low pressures. 

The glass beads were cleaned and dried to avoid any 
agglomeration (electrostatic forces or moisture). The 
glass beads were then deposited into a flexible container 
(Tygon sleeve) of 3 cm height and 2.5 cm radius. Trans- 
ducers and a pair of linear variable differential trans- 
formers (LVDT, for measurement of displacement) were 
placed at the top and bottom of the flexible membrane 
(see Fig. 0. 

Before starting the measurements, a series of tapping 
and vibrations were applied to the container in order to 
let the grains settle into the densest possible packing. 
Our goal is to establish the sample in the reversible state 
described by, e.g. Figure 2 of Nowak, et al 39]- The 
entire system was then put into a pressure vessel filled 
with oil (see Fig. . We then applied confining pressures 
ranging from to 140 MPa. The pressure was cyclically 
applied several times until the system exhibited minimal 
hysteresis. At this point shear and compressional waves 
were propagated by applying pulses. The sound speeds 
and corresponding moduli were obtained by measuring 
the arrival time from "head to head" of the transducers 
for the two sound wave types. 



Acoustic measurements 



The results we obtain are plotted in Fig. [2 and they 
compare well with the available data of Domenico [21| 
for the range 0-40 MPa. There remains an hysteresis 
component between the cycle upwards and downwards 
in pressure which is representative of the packed system. 
A more detailed comparison with theory and simulations 
is done in Figs. (0J and JSJ, below. 

Because of the deformation of the glass beads the 
height of the system decreases with the increasing pres- 
sure. In order to obtain the correct velocities from the 
arrival time of the signal, we accurately measure the dis- 
placement of the transducers as the pressure is increased 
with a pair of LVDTs. In order to avoid fracture of the 
particles due to the external pressure we use small parti- 
cle sizes to reduce the intensity of the contact forces. To 
get a qualitative idea of the presence of crushing within 
the system we observe the sample under a microscope 
after the experiments. We find that crushing occurs only 
in a very small fraction of the beads. When the exper- 
iment was repeated with larger beads of diameter 0.3 
mm many particles appeared to be crushed after apply- 
ing pressures of 140 MPa. Moreover, during this test 
there was a strong acoustic emission and a severe inflec- 
tion in the sound speeds could be noticed as we increased 
the pressure during the first cycle upwards due to the 
crushing of the beads. For the beads of size 45/im, no 
inflection is observed and no acoustic emission is heard 
during the experiment. 



IV. NUMERICAL SIMULATIONS 

We perform MD simulations of a system of 10000 
spherical particles in a periodically repeated cubic cell 
of approximately 4mm sides. The particles interact via 
Hertz-Mindlin contact forces and we choose typical val- 
ues for glass beads for fi g = 29 GPa and v g = 0.2 for a 
close comparison with experiments. We assume a distri- 
bution of grain radii in which Ri = 0.105 mm for half the 
grains and i?2 = 0.095 mm for the other half. Our results 
are quite insensitive to the choice of the size distribution. 
We include viscous damping terms to allow the system 
to relax toward static equilibrium as discussed in Section 

ITTeI 

The general scheme of the simulations is as follows: 
The simulations begin with a gas of 10000 grains dis- 
tributed at random positions inside the cubic cell. We 
first apply a compression protocol so that a dense ran- 
dom packing is generated corresponding to a predeter- 
mined value of the pressure. Then, an incremental in- 
finitesimal compression or shear is applied to the unit 
cell and the change in stress is computed, once the sys- 
tem re-equilibrates. Thus, we obtain the bulk and shear 
moduli for the system at each confining pressure. 



A. The reference state: numerical protocol 

One of the critical issues in this study is how to ob- 
tain a proper rigid frame of reference Eq. 0, {R}, from 
where we could calculate the elastic moduli. Our cal- 
culations begin with a numerical protocol designed to 
mimic the experimental procedure used to prepare dense 
packed granular materials at a given confining pressure. 
In the experiments, the initial bead pack is subjected to 
mechanical tapping and ultrasonic vibration in order to 
increase the solid phase volume fraction, as discussed in 
the previous section. 

During the numerical preparation stages we turn off 
the transverse force between the grains (fc t = 0); because 
there are no transverse forces, the grains slip without re- 
sistance and the system reaches the high volume fractions 
found experimentally during the initial compression pro- 
cess. We found that by preparing the system with fric- 
tional and elastic tangential forces, the system reaches 
states of lower volume fraction. A more complete study 
of this effect will be presented in Part II |4fJ. In the 
following calculation we concentrate on the preparation 
without friction, so that we can obtain the most compact 
states possible, mimicking our experimental procedure. 
We then restore the tangential Mindlin force and friction 
when we calculate the elastic constants. 

Starting with a set of non-contacting particles, we first 
apply a slow compression to bring the particles closer 
until a specified value of the pressure and coordination 
number is attained. This initial compression is speci- 
fied by the dashed lines in Fig. If the compression is 
stopped just before reaching a volume fraction of random 
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FIG. 3: Coordination number versus pressure obtained in the 
simulations, (a) Frictionless packs in 3D. The system becomes 
isostatic as p — * and Z ~ 6. (b) Frictional packs in 3D. Is 
the isostatic limit Z c = 4 reached asymptotically as p —* 0? 
See Part II |4(i|j for details, (c) Frictionless packs in 2D. Here 
the system is isostatic with Z c ~ 4 as p — » 0. 
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close packing (specified as Point A in Fig. and the 
system is allowed to relax, then system will relax to zero 
pressure and zero coordination number, since it cannot 
equilibrate below the maximum close packing fraction. 
This is indicated in Fig. [jji as the decrease of the coor- 
dination number and pressure towards zero. The com- 
pression is then continued to a point above the critical 
packing fraction at a target pressure, pt- The target pres- 
sure, is maintained with a "servo" mechanism 37] which 
constantly adjusts the applied strain rate e until the sys- 
tem reaches equilibrium at p t according to the following 
prescription: 



g(p-pt), 



(28) 



where p is the actual pressure of the system and g is a 
gain factor which is tuned to achieve equilibrium at every 
given pressure in an optimal way. 



B. Coordination number 

The above protocol is repeated for different target pres- 
sures and we obtain the average coordination number Z 
of these equilibrium states as a function of the pressure, 
as seen in Fig. Et. Several important points can be seen 
from this plot. Firstly, the average coordination number 
increases with the pressure as expected. Secondly, we 
find that the coordination number of the pack approaches 
a critical minimal value close to Z c « 6 as p — > 0. At low 
pressures, compared to the shear modulus of the beads 
(p <C 26GPa), the system behaves more like a pack of 
rigid balls. At this point the beads are minimally con- 
nected at Z c « 6, while in two dimensions (see Appendix 
E) the same preparation protocol gives Z c ss 4 (Fig. |3t). 

Such low coordination numbers can be understood in 
terms of simple constraint arguments for a system of N 
frictionless rigid particles in D dimensions El Hi El- 
We need to determine ZN/2 normal forces with DN 
equations of force balance. We find a critical coordi- 
nation number for which the equations of force balance 
are soluble as Z c = 2D. For large values of the confining 
pressure more grains are brought into contact, and the 
coordination number increases from its minimal value re- 
quired for stability; the system is underconstrained. Em- 
pirically, we find in 3D 



Z{p) = Z c 



P 



10 MPa 



0.30(5) 



(29) 



The pressure 10 MPa is significant since it determines 
the characteristic pressure of the crossover from the min- 
imal coordination number to a larger one. 

We also measure the volume fraction as a function of 
pressure and find that it approaches a critical value of 
4> c ~ 0.63 in the rigid ball limit as p — > 0: 



0(P) 



( P ) 

V 14 GPa/ 



0.62(6) 



(30) 



The value of <j> c = 0.63 « </>rcp corresponds to the vol- 
ume fraction at random close packing (RCP)jthe dens- 
est possible random packing of hard spheres |4ll II2I liflj , 
since the hard sphere limit in our system of deformable 
particles is achieved when the pressure (deformation) 
vanishes (RCP is only achieved asymptotically in our 
simulations). The exponent 0.62 is consistent with di- 
mensional arguments which would predict a value inverse 
of the power law between the force and displacement in 
the Hertz law, i.e. a 2/3 exponent. The exponent in Eq. 
(|29|l is determined by the behavior of the pair distribu- 
tion function near jamming. These exponents agree with 
similar calculations done by O'Hern et al. J44) , and they 
will de discussed in more detail in Ref. |4Cj . 

The low value of Z c is very significant (this number 
should be compared, for instance, to Z = 12 for a FCC 
packing) because at this minimal coordination the equa- 
tions for the force distribution can be solved without ref- 
erence to the state of strain in the system. This is the 
isostatic limit and the starting point of recent the- 

ories of stress distributions in granular packs 0, 0, EH • 
Concepts such as fragility and marginal rigidity depend 
on the existence of this minimally connected state. In the 
conclusions we will come back to discuss this issue. As 
previously reported in El an d H^, Eq. I|29|) provides 
a numerical evidence of the existence of the minimally 
connected state in frictionless granular packs. For other 
numerical work see |47j . 

To test the robustness of these results, we have em- 
ployed a second protocol in which the system is prepared 
by compressing to a point beyond the RCP fraction, then 
letting the grains relax to equilibrium without the servo 
mechanism. The final Z(p) curve is essentially identical 
to the one shown in Fig. For this reason we be- 

lieve that we have accurately approximated the reversible 
state of dense random packing, in the sense discussed by 
Nowak, et al. [39]]. 

It is important to recall that the above results have 
been obtained for a system without friction. A similar 
preparation protocol for grains with friction gives rise 
to different packings with lower coordination number. 
Similar constraints arguments as explained above give 
Z c = D + 1 for this case. Fig. |2t> shows Z(p) obtained 
for a system with friction showing that a minimal Z c « 4 
in 3D may be approached asymptotically as p — > 0, al- 
though at a slower rate than in the frictionless case. Is 
the Z c = 4 isostatic limit achieved as p — > 0? We have 
given a positive answer to this question in |46l ]. How- 
ever, recent studies 0] suggest that this may not be the 
case. We refer the interested reader to Part II of this 
work 40] for our more recent results showing that the 
rate of compression (analogous to the rate of cooling of 
a a glass-forming liquid below the glass transition) plays 
a significance role in achieving the isostatic limit in fric- 
tional packs. From now on we will concentrate on the 
calculation of the elastic properties of granular media us- 
ing the states depicted in Fig. as our starting point. 
Of course, we will restore fc t ^ for the calculation of 
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the moduli. 



C. Calculation of elastic moduli with MD 

Consider the calculation of the elastic moduli of the 
system as a function of pressure. Beginning with the 
equilibrium states of Fig. |^ we first restore the trans- 
verse component of the contact force by setting k t ^ 0. 
We then apply a small perturbation to the system and 
measure the resulting response. We don't expect slip- 
page to occur since we apply infinitesimal strain pertur- 
bations, but since we deal with a finite system we set 
the friction coefficient /jj to a large value to avoid slid- 
ing at the contacts. The elastic moduli are calculated 
by applying a given affine infinitesimal strain perturba- 
tion Ae as given by Eq. (|10fl and then monitoring the 
response of the corresponding stress a(t) as a function 
of time. After the system equilibrates again as t — > oo, 
the moduli are obtained from Eqs. (|22[) and (|24J) as the 
change in stress between the final state and the stress 
before the perturbation Aer/Ae. The procedure is re- 
peated for Ae — > to guarantee that we are testing the 
linear response regime where the elastic moduli become 
independent of Ae. Interestingly we find that the region 
where the elastic constants are well defined decreases as 
the pressure decreases. This is in agreement with the 
prediction of the EMT for the 3rd order elastic constants 



which arc found to diverge as ~ e 



-1/2 



P 



-1/3 



The shear modulus is calculated from a simple shear 
test (Aei2 = Ae 2 i ^ 0) as given by Eq. Q) 



1 Aa 12 



2 Ae 12 ' 

and also from a pure shear test with Aen = — Ae 2 2 



(31) 



1 (Ao- 22 - Aan) 

2 (Ae 22 - Aen) ' 



(32) 



We find that the values of /i determined from these 
two methods agree with each other, as expected for an 
isotropic system. 

The bulk modulus is obtained from a uniaxial com- 
pression test along the 1-direction and keeping the strain 
constant in the other directions Ae 22 = Ae 33 = 0, and 
Ae n + 0: 



K 



Ao-ii 
Aen 



(33) 



Here the stress, <7y, is determined from the measured 
forces on the grains Eq. l|2*7|l. and the strain, e^, is de- 
termined from the imposed dimensions of the unit cell. 
For instance e n = AL/L where AL is the infinitesi- 
mal change in the 11 direction and Lq is the size of the 
reference state at the given p. 



The results of our numerical calculations for K(p) and 
fj,(p) are shown in Fig. g] These results have been ob- 
tained for packings of 10000 particles. Calculations done 
with 432 spheres show similar values indicating that the 
results are free of finite size effects. We see that our ex- 
perimental and numerical results are in reasonably good 
agreement. Also shown are data measured by Domenico 
[2l|. Clearly, the experimental data are somewhat scat- 
tered at low pressure. It reflects the difficulty of the 
measurements, especially at the lowest pressures where 
there is a significant signal loss. Nevertheless, our cal- 
culated results pass through the collection of available 
data. It should be noted that the experiments are com- 
pared against the numerical results without resorting to 
the use of fitting parameters, since all the constants char- 
acterizing the grain material (/i g and v g ) are known from 
the properties of the grains. 



D. Breakdown of the EMT. Problems with fi 

Also shown in Fig. glare the EMT predictions Eq. jT2jl 
and (|14fl using the same parameters as in the simulations. 
We set Z = 6 and <f> = 0.64, independent of pressure. At 
low pressures we see that K is well described by EMT. 
At larger pressures, however, the experimental and nu- 
merical values of K grow faster than the p 1 ^ 3 law. The 
situation with the shear modulus is even less satisfactory. 
EMT overestimates flip) at low pressures but, again, un- 
derestimates the increase in flip) with pressure. 

To investigate the failure of EMT in predicting the 
correct pressure dependence of the moduli, we re-plot 
the moduli divided by p 1 ^ 3 in Fig. \5\ 

For such a plot, EMT predicts a horizontal straight line 
but we see that the numerical and experimental results 
are clearly increasing with p. It is tempting to try to 
fit the data with another power law. However, we must 
first include the power law dependence of the coordina- 
tion number and the volume fraction with the pressure 
as given by Eqs. I|29|) and l|30|) . Thus we modify Eqs. 
(|T2)l and ((Hfl to include the pressure dependence Z(p) 
and 4>{p) (this latter is a much smaller effect, see below). 
The corrected EMT is also plotted in Fig. [5] and we see 
that it predicts the same trend with pressure as the simu- 
lations. The experimental data also seem to be following 
this trend but more data over a larger pressure range are 
clearly needed. 

When analyzing K(p), we find that the corrected EMT 
is in essentially exact agreement with our numerical simu- 
lations and experimental data. Thus, we tend to conclude 
that the anomalous scaling found in the experiments is 
be a measurement of a crossover behavior as obtained by 
combining Eqs. I|12|) and i|14fl with the nonlinearity of 
Eqs. I|29l) giving rise to two distinct scaling regimes: 



K{p) ~ flip) ~ p 1/3 , for p < 10 MPa, (34) 
Kip) ~ flip) ~ p 5/9 , for 10 MPa < p < 14 GPa. 
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FIG. 4: Pressure dependence of the elastic moduli, (a) bulk 
and (b) shear moduli from MD, our experiments, Domenico 
experiments and EMT. 



Here we have not included the pressure dependence of 
the volume fraction Eq. (|30(l since it appears at the very 
large pressures above 14 GPa. At these pressures the 
beads are not supposed to follow anymore the Hertz law 
(and they may, in fact, fracture). Therefore we exclude 
this regime from our scaling analysis in Eq. (|34|l . 

Since the experiments are usually done near the 
crossover pressure of 10 MPa, it holds to reason that 
they could be measuring a crossover behavior rather 
than a true scaling regime. Moreover, even for pres- 
sures larger than 10 MPa the Hertz contact mechanics 
approach might fail since the Hertz theory is based on 
small perturbations. Thus, the true final scaling regime 
Eq. (|34|l might not be accessible experimentally, at least 
for glass beads and other rigid materials. It would be 
interesting to see if such a crossover could be observed in 
softer materials. 
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FIG. 5: Elastic moduli, (a) bulk and (b) shear, normalized 
to p 1 ' 3 and corrected EMT taking into account the pressure 
dependence of Z(p) from Fig. as well as 4>(p). 



The substitution of Eqs. an d CHH into 03 and 

(|14fl is something of an ad hoc procedure; Eqs. (|12fl and 
(|14fl were derived under the assumption that Z and <p 
are stress-independent quantities. Within the context of 
the affine assumption, the EMT derivation can be modi- 
fied to account for a continuously changing coordination 
number, Z(p). Let us assume that, in the limit of zero 
pressure, there is a probability distribution P(h) of gap 
sizes, h, between each ball and its neighbors: 



P(h) = Z c S(h) + ax + a 2 h + ... (35) 

where Z c ~ 6 represents the coordination number at zero 
stress and the rest is a Taylor's series expansion around 
h = 0. It is straightforward to re-do the derivations lead- 
ing to Eqs. I|12fl an d (|14|l following the prescription in, 
e -g- US- The results, expressed in terms of the static 
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FIG. 6: Ratio Kj \i for MD, experiments, and EMT (pressure 
independent) 




FIG. 7: K(a) and p{a) versus a for a fixed p — 100 KPa as 
calculated numerically with MD (noted MD), as calculated 
numerically using only the affine motion (noted MD AM) and 
as predicted by the EMT (noted EMT). 



compressive strain, e < 0, are 



6ir 



K = 



12tt 



Z c {-e)^ + ^ ai R){-e)^ 



(36) 



(37) 



Using a judiciously chosen value of ai ^ 0, and neglecting 
02 and all higher-order terms, a cross-plot of Eq. I|37() 
against (|36J) mimics the molecular dynamic simulations 
in Figure 0] We note, however, that, taken literally, Eq. 
(|35|l predicts Z — Z c oc p 2 ^ 3 for small p, in contrast to 
Eq. (H3)- 

Since the bulk modulus is approximately described by 
the corrected EMT, throughout the rest of the paper we 
focus on n(p). In Fig. [5]it is shown that even though the 
pressure trend is well described by the corrected EMT, 
the theory still overestimates the value of the shear mod- 
ulus. We will see later that the overestimation depicted 
in Fig. [5] becomes enormous when the tangential forces 
are diminished towards zero. In this limit, the breakdown 
of the EMT is clearly established. 

Another way of seeing the breakdown of EMT is to 
focus on the ratio K/fi, which is independent of pres- 
sure in the theory Eq. i|17|) . the simulations, and ap- 
proximately so in the experiments, as seen in Fig. [5] 
[The variation at low pressure may reflect the difficulty 
in propagating sound at low confining pressures]. The 
experiments give K//i i=s 1.1 — 1.3. Our simulations give 
K/ /i ss 1.05 ± 0.05 in good agreement with experiments. 
Notice, however, that the EMT predicts K/fi — 0.71, as 
mentioned earlier. Moreover, the effective Poisson ratio 
from simulations, v sa 0.27 is in excellent agreement with 
that of the experiment v w 0.28, but greatly differs from 
the theoretical prediction v e = 0.02, Eq. (|18fl . 



E. Role of transverse forces and rotations 

To understand why /i is overestimated by EMT we 
must examine the role of transverse forces and rotations 
in the relaxation process of the grains. These effects do 
not play any role in the calculation of the bulk modu- 
lus. According to the EMT, the transverse force F t con- 
tributes only to the shear modulus and not to the bulk 
modulus [see Eqs. (|P2)) . ljT3|l and We are there- 

fore motivated to examine the behavior of the moduli as 
a function of the strength of the transverse force. We 
replace the tangential stiffness k t in Eq. (0} by akt and 
redefine the transverse force as 



AF t =ak t (R0 1/2 As, 



(38) 



a = is appropriate for frictionlcss coupling (perfect 
slip), whereas a = 1 describes the fully frictional result 
(perfect stick) and corresponds to the results described 
so far. To quantify the role of the transverse force on 
the elastic moduli, we calculate K(a) and n(a) varying 
a from to 1, at a given pressure, p — 100 KPa, low 
enough so that the changing number of contacts does 
not play a role. 

The results arc plotted in Fig. [7| (curves labelled MD). 
To compare with the theory we plot the prediction of the 
EMT Eqs. (^2) an d 031 m which k t is rescaled by ak t 
(curves labelled EMT). (The curves labelled MD AM are 
discussed in the next subsection.) The simulation con- 
firms that K is essentially independent of the strength of 
the tangential force; both theory and simulations show a 
flat line in Fig. [7] Surprisingly, the shear modulus is ex- 
tremely sensitive to the tangential force and becomes neg- 
ligible small in the limit of frictionlcss particles (a — ► 0) 
dropping to less than 10% of the predicted EMT value. 
We see that the EMT badly fails in accounting for the 
vanishing of the shear modulus as a — > 0. By contrast 
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the bulk modulus agrees reasonably well with EMT re- 
gardless of whether there was perfect slip or perfect stick. 
What is the most serious problem with the elastic the- 
ory? In the next section we will focus on the role of 
stress relaxation and the nonaffine motion of grains due 
to disorder. 

First, however, we wish to eliminate a conceptually 
simpler effect of disorder as the explanation for the be- 
havior of fi(p). In the simulations (and presumably in 
the experiments), it is not true that each grain has the 
same number of contacts. Rather, there is a distribution 
of contacts ranging from Z = 3 to Z = 10 with a peak 
at Z = 6, which is near the average (Z = 6.14 at 100 
KPa). Thus, the local elasticity moduli can vary widely 
from one grain to another. There is a well-developed 
theory for just such situations |48|, which is also called 
a "self-consistent effective medium approximation" (sc- 
ema). Let Ki and [ii be the moduli for spherical inclu- 
sions whose volume fraction is q. The effective elastic 
constants for the composite, K* and fi*, are determined 
by the simultaneous solution of the following coupled 
equations: 



E 



K* — Kj 
Ki + (4/3)/i* 



= 0, 



and 



where 



Mi 



\ - H 



= 0, 



F* 



lx*(9K* + 8/J,*) 
6(K* + 2(i*) 



(39) 



(40) 



(41) 



Effective medium theories of this sort generally work well 
in situations in which the disorder is not too great (such 
as when there is a log-normal distribution of constituent 
properties, or when one is near a percolation threshold). 
Moreover, the sc-ema have certain desired properties, 
such as correct limitin g v alues and lying within upper 
and lower bounds. See 48] for details. 

Here we take the view that the system is a compos- 
ite consisting of spherical inclusions, each of which has 
moduli given by Eqs. I|12l) and i|14fl . In the case at hand 
it is useful to rewrite them in terms of the local value of 
the compressive strain, €{ < 0, within each inclusion [See 
Ref. [13 for details]: 



K, 



12tt 



Zi{-ei)i 



<f>(k„ + \akt) 
2(hr 



Zi(-€i)i 



(42) 



(43) 



Of course, grains with a large number of contacts, Z^, 
can be expected to have a smaller than average compres- 
sive strain, e^. In order to relate ti to the macroscopic 



strain, e*, we recognize that the spirit of the sc-ema is 
that each spherical inclusion is surrounded by the host 
material. Therefore, it is a simple elasticity problem 
to show that the differential change in strain within a 
spherical inclusion is related to the differential change in 
macroscopic strain by 



K* + (4/3)/i* 
Ki + (4/3) M * 



de* 



(44) 



We take the distribution of contacts {ci} from our sim- 
ulation at 100 KPa. It is straightforward to solve the 
system of equations (j2HJ>-(G2j!- The EMT we have been 
discussing corresponds to Zi — ►< Z > and — > e*; for 
the case a = 0, and using the same material parameters 
as before, it may be written as 



K, 



16.2(-e*) 3 / 2 , 



He = 9.7(-e s 



\3/2 



(45) 



(46) 



where the moduli are expressed in GPa. If, though, the 
full distribution of contact numbers is used in the fore- 
going analysis the results are 



K" 



15.8(-e*) 3/2 , 



fi* =9.5(-e*) 3 / 2 



(47) 



(48) 



The point of this exercise is to demonstrate that, al- 
though the packing is obviously disordered, the effect of 
the disorder alone is quite negligible as far as the macro- 
scopic elastic moduli are concerned. Similar results hold 
for a = 1. Each grain sees, more-or-less, the same av- 
erage environment as any other. In the next section 
we investigate the effects of disorder induced relaxation, 
which, we believe is the underlying effect behind the small 
values of fi(p) we are observing. 



F. Role of relaxation and disorder 

In the EMT, we saw that if an afhne perturbation of 
the form (|10fl is applied to the system, the grains are al- 
ways at equilibrium due to the assumption of isotropic 
distribution of contacts and further relaxation of the 
grain is not significant. The response is then purely elas- 
tic. 

On the contrary, in the MD simulations (and in the ex- 
periments) after the application of an affine perturbation 
via the motion of the boundaries and grains, the beads in 
the immediate neighborhood of each grain move around, 
relative to the center grain, in a way which gives rise to 
a stress relaxation associated with these rearrangements 
of particles. 
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FIG. 8: Relaxation of the shear stress (B— >C) after an affine 
motion (A— »B) in the calculation of the shear modulus. 
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FIG. 9: Shear modulus versus pressure for frictional (a — 1) 
and frictionless (a = 0) particles. 



Figure |H1 shows the behavior of ai2(t) = e\2G{t) as per 
Eq. i|20|) for a system at p — 100 KPa and with a = 
0.2 during and after the application of the affine strain 
perturbation Aei2 which moves all the grains according 
to the external strain Eq. (|10fl . We see how the system 
behaves as a viscoelastic solid as explained in Section 
III Dl When the affine perturbation is applied, the shear 
stress increases (from A to B in Fig. [SJ and the grains 
are far from equilibrium since the system is disordered. 
This is the instantaneous elastic response. The grains 
then relax towards equilibrium as (from B to C), and we 
measure the resulting change in stress Acti 2 as t — > oo 
from which the modulus /i is calculated as in Eq. (|22[l . 

For a better understanding of the approximations in- 
volved in the EMT, suppose we repeat the MD cal- 
culations now taking into account only the affine mo- 
tion of the grains and ignoring the subsequent relax- 
ation. The resulting values of the moduli are obtained 
as /i aEnc = Acrff in 7Aei 2 with Aaff nc defined in Fig. 
In Fig. 0we plot the moduli calculated in this way as a 
function of a for p = 100 KPa (curves labelled MD AM). 
The affine moduli are very close to the EMT predictions: 
there remains a 10% difference between the EMT and 
the MD (affine) which is representative of the disordered 
packing which is averaged in the EMT. Thus, the dif- 
ference between the MD and EMT results for the shear 
modulus lies mostly in the non-affine relaxation of the 
grains; this difference is largest when there is no trans- 
verse force. 

By contrast, grain relaxation after an applied compres- 
sional affine perturbation is not particularly significant 
and the EMT predictions for the bulk modulus are quite 
accurate as seen in Fig. 



G. The isostatic limit as a critical point 



the only variables with the dimension of pressure in this 
limit. A scaling argument would lead to 

I* ~ ^(p/fcn)". (49) 

The Hertz theory predicts r\ = 1/3, a result which we 
find to be valid at low pressure for frictional grains. In- 
deed, quite generally if one assumes that each grain-grain 
force scales as Eqs. @ and Q and if one assumes the 
arrangement of the grains, however disordered that may 
be, does not change with pressure then both moduli scale 
as in Eq. I|49|l with rj = 1/3. This argument specifically 
presupposes that e.g. the average coordination number 
does not change with pressure. 

Since there are no other constants that could reduce 
the value of ji for a — > we are lead to believe that a 
new exponent rj should describe the shear modulus for 
frictionless packs. This is an effect which lies outside the 
standard assumptions of elasticity theory, as indicated 
above. Since p < k n , then rj > 1/3. To give validity 
to our hypothesis, we plot in Fig. fi(p) for a = 1 
and a = 0. We see that a better fit to the low pressure 
behavior of (J>(p) for a = is achieved with rj = 2/3. 
Notice that we deliberately try to fit the data at low 
pressure to avoid the issue of the increasing coordination 
number. 

How can we explain the 2/3 scaling behavior? A pos- 
sible answer could be provided by a recent conjecture by 
Alexander jlaj who proposed the following scaling 

ti^k n A{jp)(jp/k n f. (50) 

where the function A(p) is determined by the geometry 
of the reference frame of rigidity, Eq. which is deter- 
mined, in turn, by the pressure. Assuming that the limit 
p — > is indeed a critical state of rigidity, then we expect 



The surprisingly small values we found for /j, as a — > 
raises several questions. We notice that k n and p are 



A(p)~ P x , 



(51) 



which would explain the anomalous scaling for the fric- 
tionless grains with A = 1/3, while for frictional grains 
we would have A = 0. 



H. Microstructure and force chains. 

The velocity of acoustic signals probes an effective 
medium which should be homogeneous at length scales 
larger than a typical correlation length of the material. 
Experimental and numerical work indicates that there 
is an internal structure at length scales ~ 10c?, where d 
is the typical size of the grains: the forces are observed 
to be localized along "force chains" carr ying most of the 
loads in the system (see Fig. [TTJ) ES Hi E3 A 
question of interest is how such a microstructure affects 
the properties of the system at macroscopic length scales 
where the elastic continuum theory is valid |2fi| . 

We want to quantify the relevance of force chains to 
the elastic moduli. We calculate the shear modulus as a 
function of a subset of forces belonging to the strongest 
forces in order to search for the backbone of grains which 
give rise to the shear rigidity of the material. Is this 
backbone determined by the force chains, or do the in- 
terstitial particles play also a relevant role to determine 
the rigidity? 

In this regard, recent calculations of Radjai et al. jS^ 
have shown that the stress ratio between shear and com- 
pression shows a "percolation-like" behavior: the forces 
larger than the average are responsible for most of the 
rigidity of the material. This was shown to be valid in 
2D. Here we follow [52| and define a ^-network which in- 
cludes only forces smaller than a cutoff force £. Then we 
redefine the stress Eq. H27[) and compute the shear stress 
only for the (^-network as 

|F|<C 

from which we obtain the shear modulus for the Q- 
network as /i(£) = Acri2(£)/Aei2 (for ( -* oo we recover 
our previous results). 

FigurclTUIshows the result of for p = 100 KPa and 
a = 1 and should be compared with Fig. 4 in jH^]. In 
contrast with the 2D results of we find no evidence 
of a bimodal distribution of forces which would give rise 
to a percolation-like behavior of the shear modulus. We 
see that the shear modulus and the coordination number 
increase continuously as we increase £. 

We also repeat the same calculations for our two- 
dimensional packings and find the same result as in 3D, 
i.e. we find no evidence of a bimodal character in the 
behavior of the shear modulus versus the force cut-off. 
The fact that we do not see the same behavior in 2D as in 
[52| might be related to the regularization scheme used in 
our MD simulations to handle the frictional forces which 
may eliminate the critical behavior found in |52j. Radjai 
et al. used a Contact Dynamics algorithm which tackle 
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FIG. 10: Behavior of the shear modulus and the coordination 
number for the (^-network. We use the packing at 100 KPa 
depicted in Fig. Illfa for this calculation. 

the non-smooth character of the interactions without any 
regularization schemes. 

Figurc lTTl shows our attempt to visualize force chains in 
3D packings (a) without friction under isotropic compres- 
sion, (b) with friction under uniaxial compression and (c) 
in 2D frictional isotropic packings. Force chains are not 
prominent in the 3D isotropic frictionless packing. More- 
over, the continuous variation of /u(C) obtained for this 
packing seems to indicate that all forces are important for 
the mechanical response to shear, and not just the larger 
forces which may be organized in force chains. However, 
force chains are prominent in the 3D packing under uni- 
axial compression and the 2D packing. 



V. THEORY: SINGLE PARTICLE 
RELAXATION 

Since the difficulty with the shear modulus is shown 
to be due to the relaxation of the particles from the 
initial uniform strain approximation, we next perform 
the simplest investigation that allows for some relax- 
ation. From the simulations, we know the rest positions 
of each of the particles, as well as the contact vectors 
dW = (xi — x 2 )/|x! — x 2 | (the vector from a particle to 
each of the particles with which it is in contact). Con- 
sider a specific particle. We make the approximation that 
when a small amplitude macroscopic strain is applied its 
contacting particles move according to the affine approxi- 
mation. The particle will experience an unbalanced force 
and an unbalanced torque. Accordingly, it will relax to 
a new position and orientation such that the net force 
and torque on it become zero. So, for the specific par- 
ticle we calculate its new position and orientation. Wc 
next calculate the energy stored within each of the con- 
tact "springs". We do this for each of the particles in 
the simulation to calculate the total stored energy due 
to the applied strain and we set this equal to the usual 
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FIG. 11: (Color online) Force chains in granular matter: (a) 
Frictional system under uniaxial compression near <f> c (from 
Percolating force chains are seen in this case. We apply 
an algorithm which looks for force chains by starting from 
a sphere at the top of the system, and following the path 
of maximum contact force at every grain. We plot only the 
paths which percolate, i.e., stress paths spanning the sample 
from the top to the bottom, (b) Frictionless isotropic system 
at p = 100 KPa in 3D. We plot only the forces larger than 
the average. Force chains seem to be tenuous and not well 
defined, (c) Force chains in a 2D frictional system. Force 
chains are clear in this case. 



expression for strain energy in order to deduce the new 
estimates for the bulk and shear moduli of the aggregate. 
This procedure is detailed below. 

Consider a particle, labeled a, which we take to be 
centered at the origin. It has z a contacts at the positions 
{d^ : q — l,z a }. Assuming that one of the contact 
points is displaced by an amount the increment in 
the intergrain force at contact q is 



pW = ^[(d^d^) • u (<?) ] + aK T [(I - d (9) d (9) ) • u (9) ], 
where Kn and Kt are given by: 



K 



N 



I - fa 



Kn 



W* 1/2 a/ 2 

2-v„ ' 



(54) 



(55) 



and £ is the normal displacement which can be related to 
the external pressure through the average affine approx- 
imation by 



3tt (1 



P 



<j>Z ptg 



2/3 



(56) 



The parameter a allows us to continuously investigate 
the crossover behavior from perfect slip (a — 0) to perfect 
stick (a = 1). 

As written, the total force on the specific particle, due 
to the sum of all the contact forces is not zero: 



(57) 



Accordingly, that particle will move to a new equilibrium 
position, X. Similarly, the net torque on the particle is 
unbalanced: 



(58) 



Accordingly, the particle will rotate through an angle 
w to a new orientation. The generalization of Eq. Ij53(l 
that takes into account the new position and orientation 
is 



F(?) = K N [(d^d^) • (u<«> - X)] + 

aK T [(I - d^d^) • (u<«) - X) - u X d^], 

(59) 

Now, the requirement that the particle is in equilib- 
rium with its contact forces, ^2 q F^) = 0, gives three lin- 
ear equations in the six unknowns, u) and X. The require- 
ment that the total torque must vanish, ^ d^ xF' ? ' = 



0, gives the remaining three. It is straightforward to solve 
these equations numerically. 

Having determined the new equilibrium position and 
orientation, one can show that the total work done by 
the contact forces on the a-th particle is simply 



W a = H**E£=i(d (9) -u«) a + 

aK T Eg=i I d (9) x U W | 2 -F„ • X - N„ • w}, 

(60) 

X and u} are determined as described above. In order 
to calculate W a we make the affine assumption, that the 
displacement at the contact point is simply related to the 
macroscopic strain by Eq. (|10|l . Since we know the exact 
positions of each contact vector, d 9 , from the simulations, 
we are able to evaluate Eq. (|qTJ)1 for each particle in the 
ensemble. 

We now evaluate ^ Q W a /V for a pure compression and 
for a simple shear numerically and we equate the result 
to the elastic energy, Eq. ©, in order to deduce the 
values of K and [i. 

The above procedure can only reduce the moduli rel- 
ative to those of the effective medium prediction. If, in 
Eq. Q60[l. we assume there is no relaxation [lj = and 
X = 0] , and if we replace the sum over contacts by an in- 
tegral over a presumed uniform distribution of contact di- 
rections, wc reproduce the effective medium theory, Eqs. 
Q and CH). 

The results of such a calculation are shown in Fig. 1121 
which is to be compared to Fig. [7| The static confining 
pressure is 100 KPa. We see that, relative to the effec- 
tive medium prediction, there is a small reduction of the 
bulk modulus, which is relatively insensitive to a. There 
is a much larger reduction of the shear modulus but the 
results of the simulations for the shear modulus give val- 
ues that are even smaller still. For a = (perfect slip) 
the simulations give [i = 8 ± 3 MPa, which is essentially 
indistinguishable from zero, whereas from Figure 1121 we 
have a value of 100 MPa. We see that relaxation effects 
at the single particle level, while significant, are by no 
means sufficient to explain the effect. In the fully fric- 
tional case of a = 1 there is a reduction relative to the 
EMT but the simulation gives a value of 200 ± 10 MPa. 
(In Fig. I|12fl we have extended the calculations into the 
unphysical range of a > 1 to emphasize that there is a 
slight change of slope, relative to the EMT.) 

We are thus lead to consider a more sophisticated the- 
ory in which we explicitly account for collective fluctu- 
ations. The next step in this direction is developed in 
|53| where we introduce fluctuations in pairs of contact- 
ing particles. This theory is developed for the frictionless 
case, where the reduction in shear modulus is most dra- 
matic and for which we can derive an analytic result using 
some fairly weak assumptions. 
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FIG. 12: First order correction to EMT allowing relaxation of 
grains from the affine motion. This figure should be compared 
with Fig. □ 

VI. SUMMARY AND OUTLOOK 

Where do we go from here? We clearly need new the- 
oretical frameworks to describe the collective relaxation 
of granular materials, especially under shear and for fric- 
tionless packs. Below we give a short review of some 
of the ideas that have been proposed recently, and how 
these theories are related to our results. 



A. Elastic versus fragile matter 

We have seen that the impossibility of defining a strain 
field which is inhomogeneous at the level of the grain is at 
the root of the problems of the elastic theory: the EMT 
approach relies on the assumption of a uniform strain 
field at all scales |ll3^. 

Interestingly, recent studies [l8L llH |45| have proposed 
theories of stress transmission in granular packs which 
describe the internal stresses without resorting to the use 
of strain variables, as in elasticity theory. These groups 
argue that cohesionless grains are in a "fragile state" of 
marginal rigidity or isostatic at a minimal coordination 
number Z c and they are only able to support certain 
loads without severe rearrangements. A novel closure 
relation between stress components — for instance, the 
fixed principal axis ansatz — and not between stress and 
strain — as in elasticity — has been proposed to solve for 
the indeterminacy in the granular system |45| . 

The correct type of closure relati on ( elastic or fragile) 
is still a question of much debate |54|. although there 
are recent experiments on the single-particle Green func- 
tion measurements suggesting that the elastic framework 
might be the correct approach at large scales [33 . l55t l56| . 

In the case of collective relaxation dynamics, our re- 
sults show that the elastic formulation is erroneous in 
describing the macroscopic shear response of granular 
materials. Moreover, we find that a very small shear 
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modulus appears for frictionless packs. This shear modu- 
lus decreases towards zero as p — > 0, as <p — ► </>rcP: an d as 
the system approaches the isostatic limit of Z — > Z c = 6. 

The vanishing of the shear modulus could be inter- 
preted as a "fragile" behavior. In the limit a — > a 
packing of nearly rigid particles responds to an external 
isotropic load with an elastic deformation and a finite K, 
since the external perturbation is compatible with the 
principal axes of the stress predetermined by the prepa- 
ration history of the sample. By contrast, such a system 
cannot support a shear load (/i — > 0) without severe par- 
ticle rearrangements. Thus the granular system supports, 
elastically, only perturbations compatible with the struc- 
ture of force chains and deform irreversibly otherwise, i.e. 
it is in a "fragile" state. 

B. Jamming and melting 

Our results show that the fragile limit is approached 
as the system gets closer to RCP limit, and that at RCP 
there is a jamming transition between a liquid-like state 
and a solid- like state with a finite modulus. The ap- 
proach to the critical point is characterized by several 
power-law exponents as in a second-order phase transi- 
tion. The vanishing of the shear modulus can be un- 
derstood as a melting of the system occurring when the 
system approaches the isostatic point. This fluid like be- 
havior has similarities with melting transitions found in 



compressed emulsions, and foams [13, 03 E3 near tne 
RCP fraction. A slow relaxation time and the increase of 
the correlation length between force chains is found near 
RCP. This behavior indicates that the physics of granu- 
lar materials might be closely related to other complex 
systems undergoing jamming as proposed recently |5S| 
such as glasses, colloids, foams, and emulsions. 



C. Conclusions 

Our MD simulations are in good agreement with the 
available experimental data on the pressure dependence 
of the elastic moduli of granular packings. They also 
serve to clarify the deficiencies of EMT. Grain relaxation 
after an infinitesimal affine strain transformation is an es- 
sential component of the shear (but not the bulk) mod- 
ulus. This relaxation is not taken into account in the 
EMT. 

Clearly, there is a need for alternative theories to de- 
scribe granular packings. Recent work on stress trans- 
mission in minimally connected networks may provide an 
alternative formulation and allow a proper description of 
the response of granular materials to external perturba- 
tions. 
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VII. APPENDIX 

Appendix A: Resistance against rolling and tan- 
gential force with microslip 

Our model of the intergrain contact is based on two 
assumptions. First, we consider the no-slip solution of 
Mindlin for the tangential force, and we consider total 
slip of the contact area only when the total tangen- 
tial force exceed /ifF n . However, in reality, the con- 
tact may slip over an annular ring of the contact area 
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for any finite value of the tangential force. A general 
study for several loading histories considering that mi- 
croslip occurs, i.e., |AF t | > pfAF n , was performed by 
Mindlin-Deresiewicz [59 | an d analyzed in more detail by 
Thornton and Randall 60] . They showed that the in- 
cremental tangential force can be obtained as: AF t = 
ek t (R^) 1 ^ 2 As ± fif(l — e)AF n , where e = 1 when mi- 
croslip does not occur (|AF t | < pfAF n ) and e takes 
different values depending on the path loading history of 
loading, unloading and reloading 60]. We have done pre- 
liminary tests using this more general solution of the tan- 
gential force, and found no significant changes in compar- 
ison with the results obtained with the no-slip solution of 
Mindlin. Therefore, we have performed our simulations 
using the simpler Mindlin contact theory. Besides, the 
EMT calculations are done using Hertz-Mindlin forces, 
so that we want to use the same interparticle laws for a 
better comparison between numerics and theory. 

Second, while rotation of spherical grains is allowed in 
the simulations, it is customary to model rotations with- 
out resistance against rolling at the contacts [37}]. Re- 
garding this approximation, it should be pointed out that 
some recent studies showed that resistance against 
rolling (modeled as an elastic spring yielding rotational 
resistance k r 8 r , where k r is the rotational stiffness, and 
9 r is the relative rotation by rolling) might be relevant 
for modeling shear bands. The relevancy of rotational re- 
sistance to static packings has not been determined yet, 
and therefore, we do not include it in our studies. It 
should be noted, however, that the simulations consider 
resistance against shear given by the elastic tangential 
force of Mindlin. 

Appendix B: Damping 

Recently it has been shown that in order to incorporate 
the dissipation law leading to inelasticity at the grain- 
grain contact consistent with the Hertz contact law, a 
nonlinear force dependency on the relative velocity of the 
grains in contact has to be incorporated into the contact 
law [13. 

This dissipative part of the normal force has been de- 
termined recently by Brilliantov et al. [62^] as: 

pdiss = 2 AknR i/*£ift£ 5 (B.l) 
3 

where A is a relaxation time that depends on the viscous 
properties of the grain material, and it can be uniquely 
determined from experimental measurements of the co- 
efficient of restitution for spherical beads [3^, |U ES| • 

In our studies, we are not interested in the way the 
system approaches the equilibrium state, but only in the 
final state which is supposed to be independent on the 
type of damping used. Thus, we use the more efficient 
global damping and linear contact damping described in 
Section III El However, for dynamical studies a damping 
term as in Eq. (|B.1|I should be considered as well. 



Appendix C: Model of interaction between 
droplets 

In the case of emulsions, interdroplet forces are not 
given in terms of bulk elasticity as in Hertz theory. In- 
stead, forces are given by the principles of interfacial me- 
chanics without considering shear forces [33l |3 HE EH • 
For small deformations with respect to the droplet sur- 
face area, the energy of the applied stress is presumed 
to be stored in the deformation of the surface. Hence, 
at the microscopic level, two spherical droplets in con- 
tact interact with a normal repulsive force F n ~ RyA. 
This is the so-called Princen model (65|, where A is the 
area of deformation, and 7 is the interfacial tension of 
the droplets, and R is the geometric mean of the radii 
of the undeformed droplets. Since the area of deforma- 
tion is proportional to overlap £, then the interdroplet 
interaction is F n ^7^. 

There have been more detailed numerical simulations 
|33j to improve on this model and allow for anharmonicity 
in the droplet response by also taking into consideration 
the number of contacts by which the droplet is confined. 
Typically these improved models lead to a force law for 
small deformations of the form F n oc A b , where A is 
the area of deformation and b is a coordination number 
dependent exponent rangin g fro m 1 (Princen model) to 
3/2 (Hertz model) (see also g^). 

Appendix D: Time step 

The time step is usually chosen much smaller than 
the collision time. However, since each contact is en- 
during, the collision time is extremely large and other 
conditions must be used. Besides, the collision time for 
Hertz spheres depends on the relative velocities of the 
particle, thus it does not defined a fixed time scale 

We choose the time step to be a fraction of the time 
that it takes for a sound wave to propagate on the grain. 
Moreover, the quasi-static approximation used to cal- 
culate the Hertz force is valid only when the relative 
velocities of the particles is smaller than the speed of 
sound in the grains 62]. Thus, the characteristic time 
is to = RyJ ' pg/ ' p g . Typically, one chooses a time in- 
terval much smaller than the characteristic time, then 
At = aR^J p g / ' p g with a < 1. Typical values for glass 
beads are: p = 2600 Kg/m 3 , p g w 29 GPa, R = 0.1 mm. 
Then At should be smaller than 10~ 8 s. Thus, in order 
to perform a simulation over one second, more than 10 s 
MD steps are needed, which is obviously a very intensive 
computation. In this case, it is customary to increase the 
density or decrease the rigidity of the particles to allow 
for a larger time step to integrate the equations of motion 
over realistic periods of time. If the shear modulus of the 
grains in decreased, then it should be checked that the 
resulting stresses are several order of magnitude smaller 
than p g , thus ensuring the condition of a nearly rigid 
system even though p g is taken smaller to obtain larger 
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FIG. 13: Bulk and shear moduli for a 2D packing normalized 
to p 1 ^ 3 , EMT and corrected EMT taking into account the 
pressure dependence of Z(p) from Fig. [3J: as well as <f)(p) (see 
Eqs. iD~2l and llD3l . 



time steps. 

Appendix E: Results in 2D 

Here we show the results for the bulk and shear mod- 
ulus as a function of the pressure for a two-dimensional 
pack of spherical particles interacting via Hertz-Mindlin 
forces. The 2D simulations are done with spherical Hertz- 
Mindlin balls constraint to move in a plane. Thus the in- 
terparticle force is that of the 3D case. Our system is not 
the same as a packing of disks in 2D since the latter has 
a different interaction law between particles. Our results 
are analogous to the three-dimensional case shown in Fig. 
[S] All the conclusions regarding the moduli obtained for 
3D are valid in this case as well. 

The scaling of the coordination number is similar to 
the 3D case: 



Z(p) = Z c 



P 



18 KPa m 



0.28(7) 



(D.2) 



with Z c w 4. 

For the volume fraction we obtain 



6{p) = 4> c + ( i ) 

V 32 MPa mJ 



0.4(1) 



(D.3) 



with a critical value of <\> c w 0.835, which is the RCP 
limit in 2D. This latter exponent is in disagreement with 
a mean field prediction based on the contact law, which 
would imply an exponent 2/3 [see discussion after Eq. 
(J20J] . However, we notice the large error bar of this result 
since we have only five data points. We refer to [40] for 
a more systematic study of this problem. 



